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Numerical Methods for the Discrete Map Z a 


Folkmar Bornemann, Alexander Its, Sheehan Olver, and Georg Wechslberger 


Abstract As a basic example in nonlinear theories of discrete complex analysis, 
we explore various numerical methods for the accurate evaluation of the discrete 
map Z" introduced by Agafonov and Bobenko. The methods are based either on a 
discrete Painleve equation or on the Riemann-Hilbert method. In the latter case, the 
underlying structure of a triangular Riemann-Hilbert problem with a non-triangular 
solution requires special care in the numerical approach. Complexity and numerical 
stability are discussed, the results are illustrated by numerical examples. 


1 Introduction 


Following the famous ideas of Thurston’s for a nonlinear theory of discrete com¬ 
plex analysis based on circle packings, Bobenko and Pinkall 0 defined a discrete 
conformal map as a complex valued function / : Z 2 C C — > C satisfying 

( fn.m fn+l,m) {fn+ljn+1 fn,m+l ) _ ^ 

//i+1 ,m+ i ) (fn.m +1 fn,m ) 

That is, the cross ratio on each elementary quadrilateral (fundamental cell) of the 
lattice Z 2 is — 1; infinitesimally, this property characterizes conformal maps among 
the smooth ones. A discrete conformal map is called an immersion if the interiors 

of adjacent elementary quadrilaterals are disjoint. 

A central problem in discrete complex analysis is to find discrete conformal 
analogues of classical holomorphic functions that are immersions; simply evolving 
just the boundary values of the classical function by (|T]i would not work 0. To 
solve this problem, Bobenko 0 suggested to augment (JTJ by another equation: 
using methods from the theory of integrable systems it can be shown that the non- 
autonomous system of constraints 


r ~ (ftl +1 .m fn.m ){fn .m fn—l,m) ^ C/«,m+l fn,m )(fn ,m fn.m— 1) 

a Jn,m — —— r \~2 m — —— r . 

\Jn+l,m Jn—\,m) \Jn,m +1 Jn,m—l) 

obtained as an integrable discretization of the differential equation 


( 2 ) 


af = xf x +yfy = zf z 
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Fig. 1 Left: red dots are the discrete Z 2 / 3 for 0 St n.m $C 19; blue circles are the asymptotics given 
by 0- Right: the Schramm circle pattern of the discrete Z 2 / 3 [courtesy of J. Richter-Gebert]. 


that would define f{z) = z a up to scaling, is compatible with (JTJ. Agafonov and 
Bobenko [CQ proved that, for 0 < a < 2, the system 0 and ([2]) of recursions, applied 
to the three initial values 


/o ,0 = 0, /i,o = l, h,i=e ian /\ (3) 

defines a unique discrete conformal map Z " m = f nm that is an immersion 0 Thm. 1], 
Moreover, they showed 0 Sect. 3] that this discrete conformal map Z" determines 
a circle pattern of Schramm type, i.e., an orthogonal circle pattern with the combina¬ 
torics of the square grid, see Fig. [T] They conjectured, recently proved by Bobenko 
and Its 0J using the Riemann-Hilbert method, that asymptotically 


7 " = r 


n + im 


l+O 


(■n 2 + f 


(4a) 


with the constant 


C(i — 


r(l-f) 


(4b) 


r ( 1 + f; 

For 0 < a ^ 1, as exemplified in Fig. [I] this asymptotics is already accurate to plotting 
accuracy for all but the very smallest values of n and m. If a —> 2, however, it requires 
increasingly larger values of n and m to become accurate. 

In this work we study the stable and accurate numerical calculation of Z"; to 
the best of our knowledge for the first time in the literature. This is an interesting 
mathematical problem in itself, but the underlying methods should be applicable to 
a large set of similar discrete integrable systems. Now, the basic difficulty is that 
the evolution of the discrete dynamical system 0 r and (|2j, starting from the initial 
values 0- is numerically highly unstable, see Fig.^ 


The support of the stencils of 0 and 0 has the form of a square and a five-point 
cross in the lattice Z 2 , that is, 


1 All numerical calculations are done in hardware arithmetic using double precision. 
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Fig. 2 Numerical discrete (0 ^ n,m ^ 49): recursing from the initial values 0 by a straight¬ 
forward application of the system 0 and 0 quickly develops numerical instabilities. The color 
cycles with the coordinate m. 


fn,m +1 fn +1 .m+1 
fn,m fn +1 


and 


fn— 1 .m fn,m fn+l,mi 
fn,m— 1 


with the latter reducing to be dimensional along the boundary of namely 


/o,m+l 

fo.m , resp. f n — 1.0 /„,0 /n+1,0, 

fo.m— 1 

if = 0 or ?7? = 0. Thus, the forward evolution can be organized as follows. If f n m is 
known for 0 ^ n,m ^ N the upper index N is increased to N + 1 according to: 

to compute the boundary values /at+i,o and fo ^+ 1 ; 
to compute the row 1 ^ m ^ N; 

to compute the column f n .N+ 1> l ^n ^N; 
to compute the diagonal value fy+i^N+i- 

It is this algorithm that gives the unstable calculation shown in Fig. [2] Alternatively, 
one could use (0! to calculate the row values /j v | i , m , 1 ^ ??? ^ N — 1, and column 
values fnjt+i, 1 ^ n ^ N — 1, up to the first sub- and superdiagonal (note that these 
calculations do not depend on order within the rows and columns). The missing 
values are then completed by using 0- However, this alternative forward evolution 
gives a result that is visually indistinguishable from Fig. 0 

As can be seen from Fig. 0 the numerical instability starts spreading from the 
diagonal elements /„ „. In fact, there is an initial exponential growth of numerical 
errors to be found in the diagonal entries, see Fig.0 Such a numerical instability of 


1. use 

2. use 

3. use 

4. use 
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Fig. 3 Numerical error of the diagonal values f n u from Fig.|2|(blue), of the x n as in 0 and computed 
by forward evolution of the discrete Painleve II equation (red), and of the corresponding invariant 
\x n \ = 1 (yellow); a = 2/3. They share the same rate of initial exponential growth. 

an evolution is the direct consequence of the instability of the underlying dynamical 
system, that is, of positive Lyapunov exponents. 



As a remedy we suggest two different approaches to calculating Z" m . In 


we stabilize the calculation of the diagonal values by solving a boundary value 
problem for an underlying discrete Painleve II equation and in Sections [3}|8] we 
explore numerical methods based on the Riemann-Hilbert method. The latter reveals 
an interesting structure (Section |4|: the Riemann-Hilbert problem has triangular 
data but a non-triangular solution; the operator equation can thus be written as a 
uniquely solvable block triangular system where the infinite-dimensional diagonal 
operators are not invertible. We discuss two different ways to prevent this particular 
structure from hurting finite-dimensional numerical schemes: a coefficient-based 
spectral method with infinite-dimensional linear algebra in Section |6]and a modified 
Nystrom method based on least squares in Section [8] 


2 Discrete Painleve II Separatrix as a Boundary Value Problem 

Since the source of the numerical instability of the direct evolution of the discrete 
dynamical system G and G is found in the diagonal elements /„ „, we first express 
the f n ,n directly in terms of a one-dimensional three-term recursion and then study 
its stable numerical evaluation. To begin with, Agafonov and Bobenko JT] Prop. 3] 
proved that the geometric quantities 



argx„ G (0,7 t/2), 


(5) 
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have invariant magnitude \x n | = 1 (see the circle packing in Fig. |Tj> and that they 
satisfy the following form of the discrete Painleve II equation 




1) 


f %n+ 1 
\ i XnXn+ 1 


-n{xl + l) 


( x„-i+ix„ 
\i+x n - \x n 


— ux n . 


( 6 ) 


with initial value xq = e mK / 4 . Note that for n = 0 this nonlinear three-term recurrence 
degenerates and gives the missing second initial value, namely 

_ *o(*o + q-i) 

1 i((a — 1),Yq + 1) 


Reversely, given the solution x n of this equation, the diagonal elements /„ „ can be 
calculated according to the simple recursion Q] p. 176] 

Mn = j r n+l = Mn ' ImX n , gn+1 = gn 7n-(-l,n-{-l = 8n+l& ^ i (8) 


with inital values go = 0, ro = 1 (note that u„, r n , g„ are all positive); the sub- and 
superdiagonal elements f n +i,n an d fn.n+\ are obtained from (JTJ) and ([5]) by 


fn +1 .n — 
fn,n +1 = 


( x n l)fn,n + (x„ + l)/n+l,n+l 

2xl 

(1 — x n)fn,n + (1 +' t »)/n+l,n+l 


(9a) 

(9b) 


However, given that x n is a separatrix solution of the discrete Painleve II equation Jl] 
p. 167], we expect that a forward evolution of (|5J, starting with the initial values xo and 
xu suffers from exactly the same instability as the calculation of the diagonal values 
fn,n by evolving 0 and (j2}. Fig. [3] shows that this is indeed the case, exhibiting the 
same initial exponential growth rate; it also shows that the deviation of the calculated 
values of \x n \ from its invariant value 1 can serve as an explicitly computable error 
indicator. 

In the continuous case of the Hastings-McLeod solution of Painleve II, which also 
constitutes a separatrix, Bomemann 0 Sect. 3.2] suggested to address such problems 
by solving an asymptotic two-point boundary value problem instead of the originally 
given evolution problem. To this end, one has to solve the connection problem first, 
that is, one has to establish the asymptotics of x n as n —> °o. By inserting the known 
asymptotics 0 of Z" m into the defining equation <[5), we obtain 


x n = e !7t / 4 (l +0(n 1 )) («—►«>). 

Since in actual numerical calculations we need accurate approximations already for 
moderately large n, we match the coefficients of an expansion in terms of n 1 to the 
discrete Painleve II equation ([6]) and get, as n —► °°, 
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Fig. 4 Numerical discrete Z n [ m (0 ^ n,m ^ 49): evolving from accurate values of f n +\m f n ,n+\ 
close to the diagonal back to the boundary by using the cross-ratio relations 0 develops numerical 
instabilities. The color cycles with the coordinate m. 


iff/ l(a—1) — o~ -\- (2— 2 z)a—(1 — 21) i (a 3 — (3 — 2 z)a~— (1 -T4z)aT- (3 -(-21)) 

A " = e4 V 16 m 5 

3a 4 — (12— 121)a 3 — (2 + 36z)a 2 + (28 + 41)a — (17 — 20/) 

+ 128n 4 

( (3a 5 - ( 15 - 12()a 4 - (30 + 48()a 3 + ( 150 + 24i)a 2 - (5 - 48z)a - ( 103 + 361)) 

+ 256 n 5 

+ 0(n- 6 ) V 


We denote the r.h.s. of this asymptotic formula, without the 0{n ~ 6 ) term, by x n $. 

Next, using Newton’s method, we solve the nonlinear system of N+ 1 equations 
in N+ 1 unknowns xq ..... x,v given by the discrete Painleve equation (|6]i for 1 ^ n ^ 
N — 1 and the two boundary conditions 

xo(xft+a— 1 ) 

*i - TV!- , 2 i n ' * *n = x n , 6 - 

i((a— 1 )xq + 1) 

Note that the value xq = e ,alt / 4 is not explicitly used and must be obtained as output 
of the Newton solve, that is, it can be used as an measure of success. We choose N 
large enough that |jrjv,61 = 1 up to machine precision (about N ~ 300 uniformly in a). 
Then, using the excellent initial guesses (for the accuracy of the asymptotics cf. the 
left panel of Fig. [T} 
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Fig. 5 Numerical discrete Z nm (0 ^ n,m ^ 49): recursing from accurate values of f n +i,n, 
fn,n+l close to the diagonal back to the boundary by using the discrete differential equation (|2j is 
perfectly stable. The color cycles with the coordinate m. 


4 0) = e mnlA , xi 0) =x, h6 (1 

Newton’s method will converge in about just 10 iterations to machine precision 
yielding a numerical solution that satisfies the invariant \x n \ = 1 also up to machine 
precision. Since the Jacobian of a nonlinear system stemming from a three-term 
recurrence is tridiagonal, each Newton step has an operation count of order 0(N). 
Hence, the overall complexity of accurately calculating the values x n , 0 ^ n ^ N, is 
of optimal order 0(N). 

Finally, having accurate values of x„ at hand, and therefore by ([8]» and (|9} also 
those of f n n , f, (. | „ and one can calculate the missing values of f n m row- 

and column-wise, starting from the second sup- and superdiagonal and evolving 
to the boundary, either by evolving the cross-ratio relations <|T]) or by evolving the 
discrete differential equation ([2]). It turns out that the first option develops numerical 
instabilities spreading from the boundary, see Fig. [4} whereas the second option is, 
for a wide range of the parameter a, numerically observed to be perfectly stable, 
see Fig. [5] Note that this stable algorithm only differs from the alternative direct 
evolution discussed in the introduction in how the values close to the diagonal, that 
is fn,n, fn+i,n and are computed. 

The total complexity of this stable numerical calculation of the array f n m with 
0 ^ n,m ^ N is of optimal order 0(N 2 ). 
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Fig. 6 Contours for the X-RHP p. 15]. Left: two non-intersecting circles J] centered at ±1 
(black); right: additional circle I 2 centered at 0 (red) for standard normalization at 00 . 


3 The Riemann-Hilbert Method 

Based on the integrability of the system 0 and (|2j, by identifying (JTJ as the com¬ 
patibility condition of a Lax pair of linear difference equations Q and by using 
isomonodromy, Bobenko and Its j4j p. 15] expressed the Z“ map in terms of the 
following Riemann-Hilbert problem (which is a slightly transformed and transposed 
version of the X-RHP by these authors): Let Ij be the oriented contour built of two 
non-intersecting circles in the complex plane centered at z = ±1 (see Fig.[6]left), the 
holomorphic function X : C\T\ —> GL(2) satisfies the jump condition 

z + (C) = g 1 (C)x 4C) (Ceil) 

with the jump matri>0 

Gl ^~ (y W2£—«/2(£ _ 1) -m (f + 1 )- 

subject to the following normalization 

( s+'t q \ 

0 z _ mr ](l + 0(z-')) (z->•«>)• (10c) 

Here, we restrict ourselves to values of n and m having the same parity such that 
(,m + n )/2 is an integer. The discrete Z" map is now given by the values f n _ m extracted 
from an LU -decomposition at z = 0, namely. 


O' 


(10a) 


(10b) 


2 To make G\ holomorphic in the vicinity of I] we place the branch-cut of £ a ' 2 at the negative 
imaginary axis, that is, we take, using the principal branch Log of the logarithm, 

_ e iajt/4 e -fLog(f/j). 
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Fig. 7 Left: Contour of the S-RHP [4j pp. 24-27] centered at z = 0. Right: Modified contour after 
normalizing the RHP at z = 0 and z —> °° (analogously to Fig.[6}; the relative size of the inner and 
outer circles is chosen depending on n and m and has a major influence on the condition number of 
the spectral collocation method. Proper choices steer the condition number into a regime, which 
corresponds to a loss of about three to eight digits, see 1241 Sect. 5.4], 



that is. 



Subsequently, using the Deift-Zhou nonlinear steepest decent method, Bobenko and 
Its J4|l transform this X-RHP to a series of Riemann-Hilbert problems that are more 
suitable for asymptotic analysis. The last one of this series before introducing a global 
parametrixr] the S-RHP (4| pp. 24-27], is based on the contour shown in the left 
part of Fig. p] This rather elaborate S-RHP is, after normalizing at z = 0 and z —> °° 
appropriately, amenable to the spectral collocation method of Olver fl8l : we skip 
the details which can be found in the thesis of the fourth author G.W. [[24] Sect. 5.4] 
that extends previous work on automatic contour deformation by Bornemann and 
Wechslberger nMa. Here, the relative size of the inner and outer circles shaping 
the contour system shown in the right part of Fig.[7]have to be carefully adjusted 
to the parameters n and m to keep the condition number at a reasonable size. The 
complexity of computing f n m for fixed n and m is then basically independent of m 
and n. 

In the rest of this work we explore to what extent the analytic transformation from 
the X-RHP to the S-RHP is a necessary preparatory step also numerically, or whether 
one can use the originally given X-RHP as the basis for numerical calculations. To 
this end, we replace the normalization (| 10c [> by the standard one, that is. 


x(z)=i+o(z ! ) (z —y °°), 


(10d) 


3 Though the parametrix leads to a near-identity RHP. the actually computation of the parametrix 
would require solving a problem that is, numerically, of similar difficulty as the S-RHP itself. 
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and introduce a further circle J2 as shown in the right part of Fig. [6] with the jump 
condition 




(10e) 


We define F = f] UI2 and put G(Q = Gj(Q for £ £ I) (j = 1,2). That way, the 
Riemann-Hilbert problem is given in the standard form 


x + (C) = G(C)x_(C) (Cer), x( z ) = /+0(z- 1 ) (z-►■»). (ii) 

Because of detG = 1, the solution X £ C ffl (C\F,GL(2)) is unique, see lfl2l p. 104], 


4 Lower Triangular Jump Matrices and Indices 

We note that the jump matrix G defined in (jTObJi and ( | 1 0e[ > is lower triangular. How¬ 
ever, even though the non-singular lower triangular matrices form a multiplicative 
group and the normalization at z —► 00 is also lower triangular, the solution X turns 
out to not be lower triangular. Arguably the most natural source of RHPs exhibiting 
this structure are connected to orthogonal polynomials. By renormalizing at z —> °° 
the standard RHP for the system of orthogonal polynomials on the unit circle with 



Though one could perform a set of transformations to this problem that are standard 
in the RHP approach to the asymptotics of orthogonal polynomials on the circle, 
basically resulting in an analogue of the S-RHP of (4), our point here is to understand 
the issues of a direct numerical approach to the X-RHP ( fTTj ) in a simple model case. 
It is straightforward to check that the unique solution of © is given explicitly by 



4 The standard form, see (2] p. 1124], of that orthogonal polynomial RHP would be 



The model problem © is obtained by putting the diagonal scaling at z —> °° into the jump matrix. 
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with 

t 2 t*-i r(k 

^>= 1 « + ii +, " + (fc lT T= e 'W' ,13) 

where F(z) and r(k,z) denote the Gamma function and the incomplete Gamma 
function. In particular, we observe that K] 2 (0) = -1 / 0. 

The nontrivial 12-component v of a Riemann-Hilbert problem with lower triangu¬ 
lar jump matrices, such as © or © can be expressed independently of the other 
components, it satisfies a homogeneous scalar Riemann-Hilbert problem of its own. 
Namely, denoting the 11-component of G by g, we get 

v+(C)=s(Ov-(C) (Cer), v (z) = o(z~ 1 ) (z-►<»). ( 14 ) 

If the contour is a cycle as in ( fTT| , or as in the model problem above, the general 
theory ifTTl §127] of Riemann-Hilbert problems with Holder continuous boundary 
regularity states that the Noether indedj K of ( [T4| is given by the winding number 

K = indj- g. 

More precisely, the nullity is the sum of the positive partial indices and the deficiency 
is the sum of the magnitudes of the negative partial indices, see E7J Eq. (127.30)]. 
Since there is just one partial index in the scalar case, the nullity of ( p~4] > is K if K > 0, 
and the deficiency is — K if K < 0. 

Thus, in the case of the RHP © the nullity of the scalar sub-RHP for the 
12-component is 


ind r s = ind r , 1 +indr 2 £ (n+m)/2 = 

in the case of the model RHP ( fl2| the corresponding nullity is m. In both cases, the 
unique non-zero solution of ( [T4| that is induced by the solution of the defining 2x2 
RHP is precisely selected by the compatibility conditions set up by the remaining 
linear relations of that RHP: the homogeneous part of these relations must then have 
Noether index — K. 


Impact on Numerical Methods 

This particular substructure of a Riemann-Hilbert problem with lower triangular 
jump matrices G is a major challenge for numerical methods. If a discretization of 
the 2x2 RHP induces a discretization of the scalar subproblem ( fl4| that results in 
a homogeneous linear system with a square matrix .Sty (that is, the same number of 
equations and unknowns), there are just two (non exclusive) options: 


5 Here, we identify a RHP with an equivalent linear operator equation Tu = ■■■, see, e.g., CD in 
the next section. We recall that A = dimkerT is called the nullity, p = dimcokerr the deficiency 
and K = A — p the Noether index of a linear operator T with closed range. 
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• Sn is non-singular, which results in a 12-component vjy = 0 that does not converge; 

• the full system is singular and therefore numerically of not much use (ill- 
conditioning and convergence issues will abound). 

Such methods compute fake lower triangular solutions, are ill-conditioned, or both. 

To understand this claim, let us denote the 12-component of the 2x2 discrete 
solution matrix by vn and the vector of the three other components by vv ; y. By inher¬ 
iting the subproblem structure such as ( [T4| for the 12-component, the discretization 
results then in a linear system of the block matrix form 

\ • T n J \w n J \»J 

=An 

Because of det(A^) = det(5jv)det(7jv) a non-singular discretization matrix Ay im¬ 
plies a non-singular Sn and, thus, a non-convergent trivial component v,y = 0. Such 
a non-convergent zero 12-component is what one gets, for example, if one applies 
the spectral collocation method of fl8l (with square contours replacing the circles) 
to the Riemann-Hilbert problems © or to the model problem ( [T4| . As a hint of 
failure, the resulting discrete system is ill-conditioned; details which can be found in 
the thesis of the fourth author G.W. l24l Sect. 5.4], 

The deeper structural reason for this problem can be seen in the fact that the 
Noether index of finite-dimensional square matrices is always zero, whereas the 
index of the infinite-dimensional subproblem ( fl4| is strictly positive. 

We suggest two approaches to deal with this problem: first, an infinite-dimensional 
discretization using sequence spaces, that is, without truncation, and using infinite¬ 
dimensional numerical linear algebra, and second, using underdetermined discretiza¬ 
tions with rectangular linear systems that are complemented by a set of explicit 
compatibility conditions. 


5 RHPs as Integral Equations with Singular Kernels 


In this section, we recall a way to express the RHP ( [TTj ), with standard normalization 
at infinity, as a particular system of singular integral equations, cf. ifrannim w e 
introduce the Cauchy transform 


C/ft) = 


r Q-z 


1 

2 ni 


(z?n 


and their directional limits C± when approaching the left or right of the oriented 
contour F, defined by 


C±/(rj) = lim — 

z-^rj± 27ll 



(rjer). 




Numerical Methods for the Discrete Map Z' 


13 


Note that C± can be extended as bounded linear operators mapping L 2 (T) (or spaces 
of Holder continuous functions) into itself, and C (suitably extended) maps such 
functions into functions that are holomorphic on C \ F, see m p. 100]. By using 
the decomposition C+ — C_ = id, the ansatz (by letting C act component-wise on the 
matrix-valued function u) 

X(z) =I + Cu(z), «ei 2 (r,C 2x2 ), (15a) 

establishes the equivalence of a RHP of the form ( fTTj ) and the system of singular 
integral equations 

(id —(G —7)C_)w = G —7 (15b) 

As the following theorem shows, singular integral operators of the form 

T g = id —(G — /)C_ : L 2 (F,C 2x2 ) —► L 2 (F,C 2x2 ) 

can be preconditioned by operators of exactly the same form. 

Theorem 1. Let r be a smooth, bounded, and non-self intersec tin ^contour system 
and G : r —¥ GL(2) a system of jump matrices which continues analytically to a 
vicinity ofT. Then, T G - i is a Fredholm regulator ofT G , that is, T G ~ i T G = id +K with 
a compact operator K : LrfT , C 2x2 ) L 2 (r, C 2x2 ) that can be represented as a 
regular integral operator. 

Proof The Sokhotski-Plemelj formula |17, Eq. (17.2)] gives that 2C- = — id +77, 
where 77 denotes a variant of the Hilbert transform (normalized as in flTl ). 

nno^.f rT7 dri 

m Jr ii — Q 

with the integral understood in the sense of principle values. This way, we have 
T G =A l id+B l H, Ai = i(/ + G), fli=*(7—G), 

r G -i =A 2 id+B 2 77, A 2 = ^(7 + G^ 1 ), B 2 =i(7-G- 1 ). 

By a product formula of Muskhelishvili [17, Eq. (130.15)], which directly follows 
from the Poincare-Betrand formula H3 Eq. (23.8)], one has 

T G - l T G =A\d+BH + K, 

where K represents a regular integral operator and the coefficient matrices A and B 
are given by the expressions 

6 Points of self intersection are allowed if certain cyclic conditions are satisfied OH: at such a point 
the product of the corresponding parts of the jump matrix should be the identity matrix. These 
conditions guarantee smoothness in the sense of 1261 . where the analog of Theorem[I]is proved for 
the general smooth Riemann-Hilbert data. 
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A=A2A\+BiBi, B=AoB\+B2Ai. 


Here, we thus obtain A = I and B = 0, which finally proves the assertion. □ 

This theorem implies that the operator To is Fredholm , that is, its nullity and 
deficiency axe finite. In fact, since in our examples detG = 1, we have that the 
Noether index of Tq is zero. The possibility to use the Fredholm theory is extremely 
important in studying RHPs: it allows one to use, when proving the solvability 
of Riemann-Hilbert problems, the “vanishing lemma” H26I , see also 112] Chap. 5]. 
For the use of Fredholm regulators in iterative methods applied to solving singular 
integral equations, see ll23l . 


6 A Well-Conditioned Spectral Method for Closed Contours 

We follow the ideas of Olver and Townsend l20l on spectral methods for differential 
equations, recently extended by Olver and Slevinsky E3 to singular integral equa¬ 
tions. First, the solution u and the data G — I of the singular integral equation ( | 1 5b[ i 
are expandetQin the Laurent bases of the circles that built up the cycle T. Next, the 
resulting linear system is solved using the framework of infinite-dimensional linear 
algebra mm, built out of the adaptive QR factorization introduced in |20l . 

To be specific, we describe the details for the model RHP ( [T2 | i, where the cycle T 
is just the unit circle. Here, we have the expansions 


«(C)=l£4C*, G(C)-i=Y, A kC k (Cer), 



both rapidly decaying with 2x2 coefficient matrices t4 and A^. In the Laurent basis, 
the operator C- acts diagonally in the simple form 



which gives 


C-u(Q = - L u -k£~ k (Cer). 


Note that —C_ acts as a projection to the subspace spanned by the basis elements 
with negative index. This way, the system ( |15b| i of singular integral equations is 
transformed tc(3 


7 It is actually implemented this way in Singular IntegralEquat ions . jl, a JULIA software 
package described in ED. 

8 We use the Iverson bracket of a condition: [SA] = 1 if the predicate & is true, [3A\ = 0 otherwise. 
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Uk+ 52 j < Q\AjUk-j — Afc (k&’L). ( 16 ) 

j =-°° 

Up to a given accuracy, we may assume that the data is given as a finite sum, 

G(C)-/« £ A*C*, 

&=—«j 

likewise for G _ 1 (£) — / with a truncation at t %2 . Thus, writing the discrete system © 
in matrix-vector form, the corresponding double-infinite matrix has a bandwidth of 
order 0{n \). Preconditioning this system, following Theorem[l] by the multiplication 
with the double-infinite matrix belonging to G 1 instead of G, results in a double¬ 
infinite matrix that has a bandwidth of order 0{n\ + 112 ). Since the right hand side 
of m is truncated at indices of magnitude Oin \), application of the adaptive QR 
factorization l20l , after re-ordering the double-infinite coefficients as Go, U -\, U \,... 
in order to be singly infinite, will result in an algorithm that has a complexity of order 
0{{ri\ +« 2 ) 2 « 3 ), where 773 is the number of coefficients needed to resolve u, dictated 
by a specified tolerance. 

Remark 1. The extension to systems F of closed contours built from several circles 
is straightforward. The jump data and the solution, restricted to a circle centered at a 
are expanded in the Laurent basis (z — a) k , k £ Z. When instead evaluated at a circle 
centered at b, a change of basis is straightforwardly computed using 

(z -ay = i Q (b a y- k ( z -b) k u e Z), 

valid for \z~b\ < \b — a\. Because of a geometric decay, one can truncate those series 
at k = G(l) as long as \z — b\ < 0\b — a\ with 0 < 0 < 1 small enough. The adaptive 
QR factorization can then be applied by interlacing the Laurent coefficients on each 
circle to obtain a singly infinite unknown vector of coefficients. 


Numerical Example 1: Model problem 

Because of the entries and C in the jump matrix of the model problem ( | 1 2| i, 
we have that 771 , 772,773 = 0(in) in order to resolve the data and the solution; hence 
the computational complexity of the method scales as 0(m 3 ). Using the JULIA 
software package SingularlntegralEquations . jlJjfvO.0. 1) the problem 
is numerically solved by the following short code showing that the user has to do 
little more than just providing the data and entering the singular integral equation 
( |15b| i as a mathematical expression: 

using ApproxFun, SingularlntegralEquations 


9 https://github.com/ApproxFun/SingularlntegralEquations.jl, cf. 03 - 





16 


Folkmar Bornemann, Alexander Its, Sheehan Olver, and Georg Wechslberger 


3 m = 100 

r = Circle (0.0,1.0) 

5 G = Fun(z -> [z A m 0; exp(z) l/z A m],F) 

6 C = Cauchy (-1) 

v @time u = (I-(G—I)*C)\(G—I) 

8 Y = z -> I + cauchy (u,z) 

9 err = norm(Y(0)-[0 -1; 1 (-1) A m*exp ( -Ifact (m))],2) 

The run timif^jis 2.8 seconds, the error of Y (0) is 4.22 • 1CP 15 (spectral norm), which 
corresponds to a loss of one digit in absolute error. 

Numerical Example 2: Riemann-Hilbert Problem for the Discrete Z 2 / 3 

Now, we apply the method to the Riemann-Hilbert problem {H} encoding the 
discrete Z“ map. Here, because of the exponents —m, —n and ±(n + m)/2 in (jTOl, 
we have ni,U 2 ,n^ = 0(n+m) in order to resolve the data and the solution, see Fig.[8j 
hence the computational complexity scales as 0((n + /n) 3 ). Note that this is far from 
optimal, using the stabilized recursion of Section [2] to compute a table including 
Z" m would give a complexity of order 0((n + m) 2 ). Once more, however, the code 
requires little more than typing the mathematical equations of the RHP. 

i using ApproxFun, SingularlntegralEquations 


3 a = 2/3 

i n = 6; m = 8; # n+m must be even 

5 pow = z -> exp(lim*a*pi/4)*exp(-a/2*log(z/lim)) 

6 r = Circle(-1.0,0.3) U Circle( + 1.0,0.3) U Circle(0.0,3.0) 

7 G = Fun(z -> in (z, T [ 3 ] ) ? [ z A ( (m+n) /2) 0; 0 1 / z A ( (n+m)/2) ] : [ 1 0; 

pow(z)/(z — 1) A m/(z + 1) A n 1],F) 
s C = Cauchy (-1) 

9 Stime u = (I-(G-I)*C)\ (G-I) 

10 X = z -> I + cauchy (u, z) ; 

11 X0 = X (0) 

ZaO = (-1) A (m+1) *X0 [2, 1] /X0 [1, 1] 

13 Zal = 3.610326860525178 + 2.568086087959661im # exact solution 

from the recursion as in Section 1.2 using bigfloats 

14 err = abs (ZaO - Zal) 

The run time is 2.7 seconds, the absolute error of Z^g is 3.38 • 10 8 , which cor¬ 
responds to a loss of about 7 digits. This loss of accuracy can be explained by 
comparing the magnitude of the 21-component of u as shown in Fig. [8] along the 
two black circles of Fig. [6] with that of the corresponding component of the solution 
matrix at z = 0, namely. 



-3.38121 -12.2073 + 8.68324/ 

12.2073 + 8.68324/ 66.0758 


10 Using a MacBook Pro with a 3.0 GHz Intel Core i7-4578U processor and 16 GB of RAM. 
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2 - 10 4 
1 ■ 10 4 
«21 0 
1 ■ 10 4 

2 l( ' ; 0. 0.5 . f.O 

t 



MO 6 



t 


Fig. 8 Left: un(Q on the circle £ = — 1 +0.3e 2lt "; right: « 2 l(?) on £ = +1 + 0.3e 2m '. The real 
parts are shown in blue, the imaginary part in yellow. Note that there are m = 6 oscillations on the 
left and n = 8 oscillations on the right; the maximum amplitude is about 1.1 • 10 4 on the left and 
7.5 • 10 s on the right. 


We observe that during the evaluation of the Cauchy transform ( |15a| >, which maps u K > 
X (0) by means of an integral, at least 5 digits must have been lost by cancellation— 
a loss, which structurally cannot be avoided for oscillatory integrands with large 
amplitudes. (Note that this is not an issue of frequency: just one oscillation with a 
large amplitude suffices to get such a severe cancellation.) 

Since the amplitudes of U 21 grow exponentially with n and m, the algorithm 
for computing Z" m based on the numerical evaluation of ( fl5j ) applied to the RHP 
GD is numerically unstable. Even though the initial step, the spectral method in 
coefficient space applied to ( | 15b[ > is perfectly stable, stability is destructed by the 
bad conditioning of the post-processing step, that is, the evaluation of the integral in 
( fT5al >. We refer to Q for an analysis that algorithms with a badly conditioned post 
processing of intermediate solutions are generally prone to numerical instability. 


7 RHPs as Integral Equations with Nonsingular Kernels 

By reversing the orientation of the two small circles in the RHP ( [TT| , and by simul¬ 
taneously replacing the jump matrix G\ by G\ = G, the RHP is transformed to 
an equivalent one with a contour system E that satisfies the following properties, 
see Fig. [9} it is a union of non-self intersecting smooth curves, that bound a domain 
Q. + to its left. By Q we will denote the (generally not connected) region which is 
the complement of Q UF. Note that the model problem ( [T2| > falls into that class of 
contours without any further transformation. 

We drop the tilde from the jump matrices and consider RHPs of the form 

<f> + (C) = G(C)<MC) (C e E), <P(z) = l + 0(z~ l ) (z-><»), (17) 
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Fig. 9 A modified, but equivalent, contour system E for the for the RHP l |l l| l obtained by reversing 
the orientation and by simultaneously replacing the jump matrix Gi by G^ 1 . Now, there is a bounded 
domain Q + (marked in green) to the left of E, cf. the original system shown in Fig. [6] 


on such contours systems E. It will either represent the aforementioned transforma¬ 
tion of E or the model problem ( [T2 | i. In particular, G is lower triangular and can 
be analytically continued to a vicinity of E. 

The classical theory developed by Plemelj (see El §126]) for such problems 
teaches the following: the directed boundary values <f> of the unique analytic solu¬ 
tion <J> : C \ E — > GL(2) of the RHP ( fTT) satisfy a Fredholm integral equation ifTTl 
Eq. (126.5)] of the second kind on E , namely 


<MO- 


1 f G-\QG{r\)-I 


2ili. 


v-C 


*-(V)drj=I (C €£) 


(18) 


understood here as an integral equation in L?(E, C 2x2 ). The matrix kernel of this 
equation, that is, 


*(C,t?) 


G- 1 (QG(t?)-/ _ c _, ( - G(t?)-G(C) 

?? - C n-C 


is smooth on E x E, since it extends as an analytic function and since the singularity 
at £ = r; is removable. Integral equations of the form ( |T8| with a smooth kernel are, 
in principle, amenable to fast quadrature based methods, see the next section. 

We note that, given the boundary values <£_(£) for £ £ E, the solution of the 
RHP (flTji can be reconstructed by 


(P(z) 


I- 

1 

. 27 li 


J_ [ 0 -(C) 
2ni Jz 

f G(C)^-(C) 




C-z 


d c 




z £ Q+- 


(19) 


In general, however, the Fredholm equation ( p~8] > is not equivalent to the RHP, see 
Ep- 387]: the Fredholm equation has but a kernel of the same dimension as the 
kernel of the associated homogeneous RHP, defined as 
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«f + (C) = g _ 1 (C)*-(C) (Cer), w( z ) = o(z~ l ) (z-►<»). (20) 


As we will show now, the kernel of the associated RHP is nontrivial in the examples 
studied in this work. 

First, we observe, by the lower triangular form of G, that the 11- and the 12- 
components of *F both satisfy a scalar RHP of the form ( fT4j i with a jump function g 
that has a winding number which is 

n + m 

md E g = - — 

for the discrete map Z a , and which is ind_r g = —m for the model problem ( | 1 2\ . Note 
that this winding number has the sign opposite to the results of Section [4] since the 
underlying 2x2 RHP is based on G~ l instead of G. Hence, the nullity of the scalar 
RHPs for the 11- and the 12-components of *F is zero and the deficiency is (n + m )/2 
(m in case of the model problem). As a consequence, the 11- and the 12-components 
of *F must both be identically zero. 

Next, since we now know that *F has a zero first row, also the 21- and 22- 
components of *F satisfy a scalar RHP of the form ( | 1 4[ ) each, but with a jump 
function g that has the positive winding number 


n + m 

indxg = — 

for the discrete map Z a , and indvg = m for the model problem ( | 1 2j i, just as dis¬ 
cussed in Section [4] Hence, the deficiency of the scalar RHPs for the 21- and the 
22-components of *F is zero and the nullity is (n + m )/2 (m in case of the model 
problem). Since both components are linearly independent from of each other, we 
have thus proven the following lemma. 


Lemma 1. The nullity of the associated homogeneous RHP ( |20| ), and hence, that of 
the Fredholm integral equation © is n + m in the case of the discrete map Z a and 
2m in the case of the model problem m 

Example 1. For the model RHP ( fl2| the smooth kernel of the Fredholm integral 
equation (p~8|) can be constructed explicitly. Here we have 


K{^t 7 ) 


G~ 1 (QG(ri)-I 

v-C 


0 


(n/cr-i 

n-C 

9 -? 9 -? 


A column of a matrix belonging to the kernel of (p~8]> satisfies the equation 



(C €-£), 


C21) 


where E is the positively oriented unit circle. We will construct solutions that extend 
analytically as u~ (z) and w_ (z) for zf=- 0, such that 
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(mz>) =™n-o«kv) (::g > ) ) 

By recalling the notation introduced in ( fl3l >) we observe, for k = 0, — 1, that 

res n=0 K n (z,ri)ri k ~ m = z k ~ m , 
res n =oK 2 i(z,n)ri k ~ m = ~z k e m ^(z), 

res v=0 K 22 (z,Ti)T] k = -z k . 

Using the coefficients a'-" 1 that induce a change of polynomial basis by 


m ~ 1 , 

Z k = Yj a kfz k e m -j (z) (k = 0,...,m-l), 
k=0 


we define the polynomials 


m— 1 


Pk (z)= L a Tj Z ] (k = 0,...,m—l), 


7=0 


each of which has degree at most m — 1. Then, the m linear independent vectors 



2r> w (0 \ 


(& = 0 ,... ,m— 1 ) 


( 22 ) 


are solutions of each. Thus, since its dimension is 2m by Lemma[I| the kernel of 
the integral equation ( fl8| is spanned by the 2x2 matrices whose columns are linear 
combinations of these vectors. □ 

The unique solution <t> of the RHP ( fTT) can be picked among the solutions 
of ( fl8] > by imposing additional linear conditions, namely n + m independent such 
conditions in the case of the discrete map Z a and 2m in the case of the model problem. 
Specifically, for the model problem ( fl2| ), we obtain such conditions as follows. First, 
since <J>_ (z) continues analytically to |z | > 1 and since <P (z) =1 + 0 (z~ 1 ) as z -A °°, 
we get by Cauchy’s formula for the Laurent coefficients at z = °° that 


1 

2ni 




[k=l\-I 


(k= 1,2,...). 


Second, by restricting this relation to the second row of the matrix <P for k = 1,..., m, 
we get the conditions 
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In fact, these conditions force all the components of the columns ( |22| > that would 
span an offset from the kernel of ( fl8| ) to be zero. 

For the Z a - RHP, similar arguments prove that the kernel of is spanned by 
matrices whose second row extends to polynomials of degree smaller than (n + m) /2 
to the outside of the outer circle in in Fig. [9] Thus, the same form of conditions as 
in ( |23| can be applied for picking the proper solution <£_(£), except that one would 
have to replace E by that outer circle and the upper index m by (n + m)/2. 


8 A Modified Nystrom Method 

Fredholm integral equations of the second kind with smooth kernels defined on a 
system of circular contours are best discretized by the classical Nystrom method m 
§12.2]. Here, one uses the composite trapezoidal rule as the underlying quadrature 
formula, that is, 



For integrands that extend analytically to a vicinity of the contour, this quadrature 
formula is spectrally accurate, see, e.g., 0 §2] or m §2]. 

Since the Fredholm integral equation ( fl8| has a positive nullity, applying the 
Nystrom method to it will yield, for N large enough, a numerically singular linear 
system. However, the theory of the last section suggests a simple modification of the 
Nystrom method: we use the conditions ( [23] ) (after approximating them by the same 
quadrature formula as for the Nystrom method) as additional equations and solve the 
resulting overdetermined linear system by the least squares method. 

Numerical Example 1: Model problem 

We apply the modified Nystrom method to the Fredholm integral equation represent¬ 
ing the model problem ( [T2| . By the sampling condition, see, e.g., 0 §2], the number 
N of quadrature points will scale as IV = 0{m ), hence the computational complexity 
scales as 0(m 3 ). To check the accuracy we compare with 



evaluated by the same quadrature formula as for the Nystrom method. For the partic¬ 
ular parameters m = 100 and N = 140 we get, within a run-time of 0.49 seconds for 
a straightforward Matlab implementation, a maximum error of these two quantities, 
measured in 2-norm, of 1.33 • 10 l4 . The condition number of the least squares matrix 
grows just moderately with nr. it is about 23 for m = 1 and about 650 for m = 1000. 
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No 


2/3 

Fig. 10 Absolute error of the approximation of Z,'„ by the modified Nystrom method vs. the 
number of quadrature points No on each of the three circles in Pig. [9] (the radii are 1/2 for the 
inner circles, 3 for the outer one); the total number of quadrature points is then N = 3 x No- One 
observes, after a threshold caused by a sampling condition, exponential (i.e., spectral) convergence 
that saturates at a level of numerical noise at an error of about 10~ 9 . Run time of a Mathematica 
implementation with No = 42 is about 0.15 seconds. 


Numerical Example 2: Discrete Z 2 / 3 

Now, we apply the modified Nystrom method to the Fredholm integral equation 
representing the RHP ( [TTj ) subject to a transformation to the form ( fT7j ). Here, the 
sampling condition requires N = 0{n +m), hence the computational complexity 
scales as 0((n + ra) 3 ). For Zgg 3 , the modified Nystrom method yields the convergence 
plot shown in Fig. [TO] it exhibits exponential (i.e., spectral) convergence until a noise 
level of about lO^is reached, which corresponds to a loss of about 6 digits. The 
reason for this loss is that this method for approximating the discrete Z" suffers the 
same issue with a bad conditioning of the post-processing step, that is, of 

as the spectral method for the singular integral equation discussed in Section [6] Here, 
the amplitude of the real and imaginary part of <£>_ (Q along the two inner circles is 
of the order 10 4 which causes a cancellation of at least 4 significant digits. 
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9 Conclusion 

To summarize, there are two fundamental options for the stable numerical evaluation 

of the discrete map Z“ m . 

• Computing all the values of the array 1 C n, m ^ IV at once by, first, computing the 
diagonal using a boundary value solve for the discrete Painleve II equation ([5]) and, 
then, by recursing from the diagonal to the boundary using the discrete differential 
equation 0. This approach has optimal complexity 0(N 2 ). 

• Computing just a single value for a given index pair ( n , m) by using the RHP © 
and one of the methods discussed in Sections [6] or [8] Since both methods suffer 
from an instability caused by a post-processing quadrature for larger values of n 
and m, one would rather mix this approach with the asymptotics 0 . For instance, 
using the numerical schemes for n, m ^10, and the asymptotics otherwise, gives a 
uniform precision of about 5 digits for a = 2/3. Higher accuracy would require the 
calculation of the next order terms of the asymptotics as in Section[2] Obviously, 
this mixed numerical-asymptotic method has optimal complexity 0(1). 
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